% 文献[2] MMSE Extension of V-BLAST based on Sorted QR Decomposition 的仿真复现
clear all;
close all;

% 仿真参数设置
Nt = 4;          % 发送天线数
Nr = 4;          % 接收天线数
SNR_dB = 0:2:30; % SNR范围
num_trials = 10000; % Monte Carlo试验次数

% 初始化BER存储
ber_mmse = zeros(size(SNR_dB));        % 传统MMSE
ber_vblast = zeros(size(SNR_dB));      % V-BLAST
ber_sorted_mmse = zeros(size(SNR_dB)); % 排序MMSE-QR
ber_mmse_vblast = zeros(size(SNR_dB)); % MMSE-VBLAST

% 添加进度显示
fprintf('开始仿真...\n');
tic;

% Monte Carlo仿真
for snr_idx = 1:length(SNR_dB)
    SNR = SNR_dB(snr_idx);
    fprintf('正在处理 SNR = %d dB (%d/%d)\n', SNR_dB(snr_idx), snr_idx, length(SNR_dB));
    
    errors_mmse = 0;
    errors_vblast = 0;
    errors_sorted_mmse = 0;
    errors_mmse_vblast = 0;
    total_bits = 0;
    
    for trial = 1:num_trials
        % 生成MIMO系统
        [y, H, x] = generate_mimo_system(Nt, Nr, SNR);
        
        % 各种检测器
        x_est_mmse = mmse_detector(y, H, SNR);
        x_est_vblast = vblast_detector(y, H, SNR);
        x_est_sorted_mmse = sorted_mmse_qr_detector(y, H, SNR);
        x_est_mmse_vblast = mmse_vblast_detector(y, H, SNR);
        
        % 计算误比特率
        [errors_mmse_trial, bits] = count_bit_errors(x, x_est_mmse);
        [errors_vblast_trial, ~] = count_bit_errors(x, x_est_vblast);
        [errors_sorted_mmse_trial, ~] = count_bit_errors(x, x_est_sorted_mmse);
        [errors_mmse_vblast_trial, ~] = count_bit_errors(x, x_est_mmse_vblast);
        
        % 累加错误
        errors_mmse = errors_mmse + errors_mmse_trial;
        errors_vblast = errors_vblast + errors_vblast_trial;
        errors_sorted_mmse = errors_sorted_mmse + errors_sorted_mmse_trial;
        errors_mmse_vblast = errors_mmse_vblast + errors_mmse_vblast_trial;
        total_bits = total_bits + bits;
    end
    
    % 计算BER
    ber_mmse(snr_idx) = errors_mmse / total_bits;
    ber_vblast(snr_idx) = errors_vblast / total_bits;
    ber_sorted_mmse(snr_idx) = errors_sorted_mmse / total_bits;
    ber_mmse_vblast(snr_idx) = errors_mmse_vblast / total_bits;
end

simulation_time = toc;
fprintf('仿真完成，总用时：%.2f秒\n', simulation_time);

% 绘制性能曲线（匹配文献格式）
figure;
semilogy(SNR_dB, ber_mmse, '-o', 'LineWidth', 2, 'DisplayName', 'MMSE');
hold on;
semilogy(SNR_dB, ber_vblast, '-s', 'LineWidth', 2, 'DisplayName', 'V-BLAST');
semilogy(SNR_dB, ber_sorted_mmse, '-^', 'LineWidth', 2, 'DisplayName', 'Sorted MMSE-QR');
semilogy(SNR_dB, ber_mmse_vblast, '-d', 'LineWidth', 2, 'DisplayName', 'MMSE-VBLAST');
grid on;
xlabel('E_b/N_0 in dB');
ylabel('BER');
legend('show', 'Location', 'southwest');
title('BER Performance Comparison (n_T = 4, n_R = 4)'); 